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1. INTRODUCTION 

Solving airfoil flow problems numerically with minimal number of 
grid points, or maximizing the accuracy of flow solutions with a fixed 
number of grid points, has become very important since more practical 
numerical solutions to three-dimensional airfoil flow problems are being 
explored. For the purpose of increasing accuracy of flow solutions, 
grid clustering techniques have been proposed [1,2]. These techniques 
were used until recently without simultaneous or iterative interaction 
with the flow solution, but rather the user arbitrarily guessed where 
fine grids might be needed and then decided how to cluster grid lines. 

As an advanced way of clustering grids, the solution-adaptive grid 
concept was lately proposed [3,4] and shown to be effective in 
increasing the accuracy of numerical solution to the transonic full 
potential equation for airfoil flow problems. In the approach of Holst 
and Brown [4],. the grid distribution along an airfoil was adapted to the 
initial flow solution. Then, grid systems in the flow field were 
obtained by the elliptic-type grid generation equations devised by 
Steger and Sorenson [2]. One drawback with this approach is that, while 
clustering grids along the airfoil is easy, maintaining the same degree 
of clustering inside the flow field is not. This is because the 
Steger-Sorenson iterative scheme becomes very slow to converge or even 
unstable when strong clustering is extended to an area far from the 
airfoil. 

This report decribes a new method of adaptive grid generation for 
the transonic full potential equation. The present approach requires 
four steps of computation: 

Step 1: Generation of initial grids 

Step 2: Solution of full potential equation on the initial grids 

Step 3: Generation of solution-adaptive grids 

Step 4: Solution of full potential equation on the adaptive grids 

*Supported under NASA-Ames University Consortium Joint Research Interchange 
NCA2-OR565-901. 
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The solution-adaptive grid generation method proposed in this 
report has some similarity to the Holst-Brown method when grids are 
adapted along the airfoil surface, but it determines grids inside the 
flow field systematically without solving the elliptic grid generation 
equation. The present approach not only saves computing time but has an 
advantage that strong clustering along the airfoil surface as well as 
inside the flow field is very easy. 

Although the present method is applicable to clustering grids in 
the 5-direction (tangent to airfoil surface) as well as in the 
n direction (normal to airfoil surface), we forcus our attention in this 
report only to that in the 5-direction. 


2. TRANSONIC FULL POTENTIAL EQUATION 

The full potential equation for air flow around an airfoil is given 

® ( 2 . 1 ) 

p = [ 1 - ( + 4>y^ )(y-1)/(y+1) (2.2) 

^bere p is the density nondimensionalized with stagnation density^ 

'j)^= u and (J) = v are velocity components, which are nondimensio- 
nalized by the critical sound speed; y is the specific heat ratio of 
air; x and y are Cartesian' coordinates. 

The boundary condition for <}) along the airfoil surface is given by 
9(j)/9n = 0 where 9/9n is the derivative normal to the surface. 

In order to apply efficiently the iterative numerical solution 
techniques such as SOR (successive-over-relaxation) or ADI (alternating- 
direction-implicit) to solving elliptic partial differential 
equations, Eq.(2.1) is transformed from the physical domain (x,y) to a 
computational domain (5,n), so the grid lines on the computational 
domain are rectangular and uniformly spaced in both directions. 

The full potential equation on the computational domain is 

(pU/J)^ + (pV/J)^ = 0 (2.3) 

where 

U = + A^^^ 

V = A^^^ + 

^ ’ ^2 “ ^x^x + ^yby . ^3 = 

^X^y ~ 

P = [1-(U<^^ + V<|)^)(y-1)/(y+1)]^^^^"^^ 
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Finite difference equations for Eq.(2.3) are written and solved 
iteratively by using the AF2 iteration scheme detailed in reference [5]. 

The TAIR (^ransonic Airf oil analysis) computer code is used for this purpose. 


3. GENERATION OF INITIAL GRIDS 

Grid generation in Step 1 (initial grids) is performed by solving 
the elliptic equations originally proposed by Thompson and et al. [1], 
and revised later by Steger and Sorenson [2]. The alternating- 
direction-implicit method has been successfully applied here to solving 
elliptic grid generation equations with grid spacing control. The 
initial grids may be generated by the hyperbolic grid generation 
equations [6], but this approach is left open for further study. 


3.1 Elliptic Grid Generation Equation and Iterative Scheme 

The grid generation equation with spacing control may be written in 
the form 


where 

A 2,2 

A ^ 

® ) 

r 2 , 2 

C = + y^ 

E = J^exp[-c (n -Ti)]Pi(0 
^ max I 

F . j2exp[-c^(n„^^-A)lQj(5) 

- V5 

In the above equations, is a user-specified constant; 
and Q^(5) are given by 


where 


= J 

Qi(5) = J 


■ ynR2>‘n= 


max 


(- y^Rj + x^R2> ln=Ti. 


max 


(3.1) 


(3.2a) 

(3.2b) 


R^ = - J 2(Ax^^ - 2Bx^^ + Cx^^)l^^^ 

-.0 max 
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The difference form of Eq.(3.1) may be written as 

Lx- + + + 

Ly - (Aa^j + 86^^ + ca^^ + Eaj + F6^)y . 0 

where 6^^ and 6^ are second order and first order difference 
operators, respectively. L in the above equations represents the entire 
difference operator after the first equality sign. 


3.2 Iterative Scheme 

The alternating-direction-implicit (ADI) Iterative scheme using 
residuals may be written as 

Mv = - (3.4a) 

Mw = - a)Ly^*^~^^ (3.4b) 

where 

V = ^ 

w = y(*') - 

and where t is the iteration number, u) is a relaxation parameter and M 
is an iterative operator that is chosen to be similar to L. As M 
becomes closer to L, the iterative scheme becomes more implicit, so the 
convergence rate increases. Here, we choose M as 

M = (a + + Fd^)(a + + Ed^)/a (3.5) 

where a is an acceleration parameter varied at each iteration. Since 
difference operators in C and n are separated, each of Eq.(3.4) can be 
solved in two steps: 

z = a (a + Adp + Edp )“^ RHS 

^ ^ (3.6) 

V or w = (a + Cd + Fd ) z 

nn h 

where RHS represents the right side of Eq.(3.4a) or Eq.(3.4b). In order 
to increase stability of iteration, either forward differencing or 
backward differencing, whichever increases the diagonal terms of the 

tridiagonal matrices, is selected for d^ and d at each grid. 

4 n ^ 

The values of P.(0 and Qi(5) are recalculated after each 
iteration by using the most updated values of x and y. However, large 
changes of P, and Q. in any two consecutive iterations cause 
instability in the iterative scheme of mesh generation. In order to 
prevent such large changes, using under relaxation for P, and Q. 
is necessary. ^ 
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4, ADAPTIVE GRID RELOCATION METHOD 
4.1 Algorithm 

Assume that an initial grid system is given, on which the initial 
solution of the full potential equation is obtained. The word 
'adaptive' indicates that a new improved grid system is generated 
adapting to the initial solution. The new grid system is obtained by 
relocating the initial grids along each ? line without changing the 
shape of the C lines. 


Let us pay attention to the grids on the airfoil surface (a 
C line). Denote s=s(C) as the circumference on the 5 line between 5 and 
5=1 (reference point) on the initial grid system (see Fig.l). The 
circumference between grid i and i=l (reference grid) may be denoted as 
Si=s(^i). Now, we move grid i to a new location on the same 
5 line at circumference Sj^'*s(5j^' ), where 5^' is the value 
of 5 for the i-th grid after moving on the Initial computational domain. 
Relocation of all the grids except the reference grid is achieved if 
one-to-one correspondence between 5 and 5' is specified by a functional 
relation as 

5 = T(5') or equivalently 5'=T ^(5) (4.1) 


In the above equation, T(5') is a monotonic function of 5', given by 


f5* 

T(5') = A[5'- 1 + <K5'')d5”] + 1 

J 1 


(4.2) 


where is the grid density function described in the next paragraph, 
and A is a normalization constant determined by the condition, 

T(5 ') = 5 =5 

max max max 

The grid density function is determined by solving 

- d^ii»(5')/d5'^ + ot(5') = f iap/35' I + glkl (4.4) 

In the above equation, a is a parameter to adjust the propagation of the 
effect of grid spacing control along the 5 line; f is the parameter that 
emphasizes the effect of the fluid density gradient: k is the curvature 
of the airfoil surface and g is to control the effect of the airfoil 
surface curvature on grid density. 

In principle, the above algorithm may be applied to each 5 line 
independently. However, we determine T(5') using p(5' ,n) along the 
airfoil surface first, and then use the same T(5') for all the 5 lines. 
The values of 5^' for new grids along the j-th n line may be denoted 
Since grids are relocated along the initial n lines, the 
valued'’of x and y for the new grids are found from the relations. 



y 


i,j 




T 
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by using iMerpolations ,^^^ w x»x('C,-n) and y«y(5yn> ar^ the s . 
transf ormatlsn between; tphe physical ahmain %nd- the initial coftputatidnal 
domain. Interpolations ate necessary, becattse xCC, n) and yC £. n) ara 
known only on the initial grids, 

4.2 Schematic Illustration 

. to explain the grid ^^elocatioh algoEitto^^ more 

qualitatively, suppose in Eq.(4.4) g-0 and Sp:/ 35 / behaves as a delta 
function at 1^, purpose, as shown in Fig. 2-A. The 

solution of Eq.(4.4) with a set of specified 0 and f is illustrated by 
a solid line in Fig. 2-B. ^ 

By introducing i|) thus determined into Eq.(4.2), T(£') with A«1 
becomes Curve-a' in Fig. 2-C. After normalization, T(£/) becomes 

Curve-a in the same figure. Curves-h and c are obtained by changing a 
and f as follows: ® ® 

Curve-b ; only f is decreased 
Curve-c : a and f are both increased 

With a decrease in f, <j> will become smaller by a factor. The 
curve-b of T(£ ) in Fig. 2-C becomes closer to a straight line which is 
the case When f=g-0. With a simultaneous increase in both a and f, the 
distribution of tj; becomes more sharply peaked about £ . 

Consequently, a very sharp increase in T(£') occurs around £ as 
illustrated by Curve-c in Fig. 2-C. ° 

Although the qualitative effect of altering a and f are explained 
by considering a delta-function-like distribution in 3p/3£', the 
distribution of 3p/3£' in actual problems is a continuous function with 
a large change across a shock. 


5. NUMERICAL STUDY 

The effect of applying the proposed adaptive grid generation is 
.tigated by considering a lifting flow of the NACA 0012 airfoil with 
the following parameters: 


Attack angle 
Number Of 


2°degree 
0.75 


grid points 99x25(coarse grid) 

The grid systems with 99x25 grid points with or without solution 
adapting will be referred to as 'coarse " 


_As the referende flow solution to which the present numerical 
results are to be compared, the same flow problem was run on a fine grid 
system (245 grids in the ^-direction and 56 grids in the q-direction) . 
The pressure coefficient distribution for the reference calculation is 
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plotted with a plain black line In every figure for the pressure 
coefficient distribution calculated with the coarse grids* The lift) 
drag and nionent coefficients calculated with the reference fine grid 
system are CL » 0.5997, CD - 0.0190, CM - -0.0268, respectively. 

Figure 3 shows the Initial grids generated in Step 1 and the CP 
(pressure coefficient) distribution obtained from Step 2. In Fig. 3-A, 
only the inner 14 C lines are shown for brevity and a better resolution 
of grid points along the airfoil surface. The line with diamonds in 
Fig. 3-B is the result of the present calculation. The diamonds show 
the location of grids. The graph is stretched toward both leading and 
^^^nlllng edges . The large disagreement in the pressure coefficient 
distribution at the shock from the reference CP distribution is typical 
when grids are coarse across a shock. 

Figure 4-A is the result of adaptive grids generated in Step 3 with 
f»10, 0-0.5. Denser grids at the leading edge are caused by large values 
of 3p/3§ at the leading edge. A strong clustering occurs also at the 
shock. In Fig. 4-B, the pressure coefficient distribution obtained 
using the adaptive grids in Fig. 4-A is compared with the reference 
pressure coefficient. The agreement in the CP distribution at the shock 
is excellent. 

Figure 5-A shows the same results as Fig. 4 except with f=10 and 
0=0.1. Clustering at the leading edge and at the shock is stronger than 
in Fig. 4-A. Figure 5-B shows that the agreement of CP distribution at 
the shock is good. However, the disagreement in CP inside the 
supersonic bubble, and after the shock, is larger than in Fig. 4-B. The 
reason is attributed to coaser grids in the 5-direction in other areas 
than the leading edge and the shock. 

Figure 6 shows the same results as Fig. 4 except with f =10 and 
0=1.5. With a larger value of o, the effects of clustering grids become 
more localized. However, when a only is increased with f unchanged as 
in this Case, the overall clustering effect becomes weaker. This is tlie 
reason why the results of Figs. 6-A and 6 -B both fall in between those 
of Fig. 3 and Fig. 4. 

In Fig. 7, clustering at the leading edge is supressed so grid 
lines only around the shock are strongly clustered. Agreement in the 
shock with the reference case is excellent. The change of CP across the 
shock is even sharper than the reference CP distribution. 

Figure 8 -A shows the adaptive grid system in which the number of 
grids along the upper airfoil surface is artificially increased while 
that in the lower surface is drastically decreased. This increase of 
grids above airfoil is done by artificially changing the value of 3 p /35 
in Eq.(4.4). In spite of finer grids in the supersonic region than the 
previous cases, agreement in the CP distribution, particularly across 
the shock, is worse. However, the shape of the CP distribution in the 
supersonic bubble is more similar to the reference case than any 
other case. The agreement in the shape is due to the fine grids in the 
5-direction, but disagreement in values is attributed to coarseness of 
grids in the n-direction. This observation suggests importance of 
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imptoving grid configuration in the leadini edge iegion of importance of 
improving the accuracy of the difference" equetiph M fhe same area# The 
disagreement of the calculated., CP diatributipii frofii the referehce 
distflbutiott in the trailing edge area is due to Coarse gflds. This tf ailing 
edge coarseness Could also cause a small effof iti the pircUiatioh 
be partly respohSible for the error in the dhoCk position 

The difference equations used fof the supersonic region are 
essentially the flrst-order-acCuate upwind differe^^^^ scheme* In order 
to maintain the upwind difference Scheme, ^Mr had aSsumed that the flow 
in the supersonic region is in the positive direction of g if §>5. 
or in the negative direction of 5 if where Ct, 

Until adaptive grids were applied j this ^gorithm did noi caMe any 
problem, because 5. was located in the vicinity of the leading edge 
any way, and maintained enough distance from the boundary of supersonic 
regions. However, when adaptive grids are USed, 5j. may move into 
a supersonic region, particularly when g fids are Strongly clustered on 
one Side of airfoil where a supersonic region appears* If this happens, 
the upwind difference scheme is violated and an instability of iteration 
is resulted* In order to Correctly apply the upwind difference scheme, 

^half value of § at the leading edge on the 

adaptive grids* Figure 9 Shows the CP distribution calculated by TAIR 

before redefiriirtg using an adaptive grid system with very strong 

clustering , where the CP distribution behaves erratically in the 
supersonic region because the upwind difference scheme was violated. 


6. CONCLUSIONS 

A new efficient method of generating adaptive grids for the 
transonic full potential equation is presented* With the present 
method, adapting grids to high curvature surfaces , shocks , and wherever 
high flow acceleration or deceleration occurs is easy. Since the 
adapting procedure needs the solution of a one- dimensional boundary 
value problem only once, the present method consumes a very small 
fraction of the computational cost required to generate a new grid 
System by the standard grid generation method* The additional advaritage 
of the present method is that an accurate initial guess for the 
iterative solution on the new adaptive grid can be easily obtained from 
the initial flow solution by Using the same interpolation procedure as 
Used for grid relocation, so computing time for the flow solution on the 
adaptive grids can be substantially saved also* 

The flow sblutidn on the adaptive grids show that accuracy in the 
pressure coefficient distribution across’ the shock is substantially 
improved with grid clustering* The CP distribution in the supersonic 
bubble is more sensitive to the grid spacing in the ^-direction than in 
the subsonic flow region* It is also affected by the grid spacing in 
the n-directioh particularly in the leading edge area where flow 
acceleration is extremely high* 

While optimal distribution of grids along the, airfoil surface is 
important, better grid clustering in the n-direction around the leading 
edge seems to be vital to further improvement of accuracy in the flow 
solution* 
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Fig. 2. A simple example using the present grid-point relocation algorithm. 
A) Density gradient forcing function. B) Resulting grid density function. 

C) Resulting clustered arc-length distribution on the airfoil surface. 
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Fig. 3-A Initial Grids 
(only 14 inner n-lines are plotted) 
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Fig. 3-B GP Distributtort 
(The graph is stretched toward 
both leading and trailing edges 
for better resioTutiohf of cuiTwes) 






f=10, 0=0.5 

Fig. 4-A Adaptive Grids 
(only 14 inner n-lines are plotted) 



Fig. 4-B CP Distribution 
(The graph is stretched toward 
both leading and trailing edges 
for better resolution of curves) 




f=iO, 0=0.1 

Fig, 5-A Adaptive Grids 
(only 14 inner n-lines are plotted) 
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Fig. 5-B CP Distributioh 
(The graph is stretched toward 
both leading and trail ihg edges 
for better resol ution of curves) 








f=10, 0=1.5 

Fig. 6-A Adaptive Grids 
(only 14 inner n-lines are plotted) 





Fig. 6-B CP Distribution 
(The graph is stretched toward 
both leading and trailing edges 
for better resolution of curves) 
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f=20, 0=0.3 

except f =0 fo r 4 0<c ' <60 


Fig. 7-A Adaptive Grids 
(only 14 inner n-lines are plotted) 
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Fig. 7-B CP Di stHbution 
(The graph i s stretched toward 
both 1 eadi ng and trai 1 i ng edges 
for better resol ution of cdrveS) 
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